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Recent work has shown that different theoretical approaches to the dynamics of the Susceptible- 
Infected-Susceptible (SIS) model for epidemics lead to qualitatively different estimates for the po- 
sition of the epidemic threshold in networks. Here we present large-scale numerical simulations of 
the SIS dynamics on various types of networks, allowing the precise determination of the effective 
threshold for systems of finite size N. We compare quantitatively the numerical thresholds with 
theoretical predictions of the heterogeneous mean-field theory and of the quenched mean-field the- 
ory. We show that the latter is in general more accurate, scaling with N with the correct exponent, 
but often failing to capture the correct prefactor. 
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I. INTRODUCTION 

Models for disease propagation are a paradigmatic ex- 
ample of processes for which the interplay of a simple 
dynamics and a topologically complex interaction pat- 
tern P, Q gives rise to nontrivial phenomena [sHs} ■ In 
this context, the Susceptible-Infected-Susccptiblc (SIS) 
model plays a paramount role, as the simplest model 
undergoing an epidemic phase transition between an ab- 
sorbing [8[, healthy phase, where infection rapidly dis- 
appears, and an active phase, with a stationary endemic 
state characterized by a finite fraction of infected individ- 
uals. In the SIS model nodes can be either susceptible 
or infected. An infected node spreads the disease to each 
one of its susceptible contacts at rate A, while it heals at 
a rate /i, usually taken to be 1. A critical value A c of the 
spreading rate separates the absorbing phase (A < A c ) 
from the endemic one (A > A c ). 

Traditional mathematical epidemiology has analyzed 
the behavior of the SIS model mainly within a mean- 
field (homogeneous mixing) approximation Q- In this 
case, the epidemic transition occurs for a finite value A c 
of the control parameter, inversely proportional to the 
average connectivity (k) of the interaction pattern. This 
view was completely revolutionized after the study of SIS 
on networks with large variations of nodes' connectivities. 
The first mathematical approaches to SIS in networks 0, 
Hol | were based on the heterogeneous mean-field (HMF) 
approach Q , an extension of standard mean- field theory 
taking explicitly into account the large fluctuations of 
the degree k (number of connections) of single nodes, and 
that is based on replacing the actual topological structure 
of the network, given by its adjacency matrix AjJ3, by 



vertices of degree ki and kj are connected in the original 
network (the so-called annealed network approximation 



From here, an expression for the epidemic threshold is 
derived for uncorrelated network^] 111 
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where (k n ) = £\ k n P(k) is the n-th moment of the net- 
work's degree distribution P(k) Eq. ((T|) predicts the 
epidemic threshold to vanish for networks with diverging 
second moment ("scale- free"), and to be finite otherwise. 
This result has huge implications, since many real net- 
works have a power-law degree distribution P{k) ~ fc -7 
with 7 < 3, so that (k 2 ) grows unboundedly with the 
system size P, Q, Q . 

For almost a decade, HMF results for the SIS model 
were considered essentially correct in the case of random 
networks within the statistical physics communitjo, al- 
though no systematic detailed investigation of their accu- 
racy was performed [14j . In other communities, however, 
different theoretical approaches have been applied, yield- 
ing opposite results. Chatterjee and Durrett [15|,[l6( have 
rigorously proven that, for strictly infinite system size, 
the epidemic threshold is exactly for any exponent 7 of 
the degree distribution. Although of fundamental impor- 
tance, this result does not provide a simple understanding 
of the physical origin of the threshold vanishing. Within 
the computer science community, another approximate 
approach, that we term quenched mean-field (QMF) the- 
ory, has been put forward by Wang et al. [l7|- The basic 
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expressing the probability that two 



A matrix taking value 1 if nodes i and j are connected and zero 
otherwise. 



2 Networks such that the probability that an edge arrives at a node 
of degree k is proportional to fcP(fc) Q{. 

3 Its validity in regular random graphs has been however contested 
by means of a HMF pair approximation theory H3 , Il3l| 
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idea is to write down the evolution equation for the prob- 
ability pi that a certain node i is infected. Taking into 
account the actual connections in the network, as given 
by the adjacency matrix, this approach predicts the exis- 
tence of a threshold given, for any graph, by the inverse 
of the largest eigenvalue of the adjacency matrix, namely 

A? MF = j-. (2) 

This result has been obtained later also by Gomez et 
al. fl8| and, in a more refined way, by Van Mieghem et 

Equation ^ can be complemented with the expression 
for Ayv obtained by Chung et al. [20( ,_to obtain an explicit 
estimate of the epidemic threshold [21| : 

A C L _ ( Ci/Vfcmax V^max > Jj?j- hi 2 (TV) 

C " I c 2 /^ 7$ > VQn(TV) ' 

where fc ma x is the largest degree in the network and C\ and 
C2 are numerical constants. This formula agrees with the 
result of Ref. [l5[ since the threshold vanishes for large 
TV as soon as fc max diverges with TV, independently of 
7 @. The threshold © scales as the HMF result only 
for 7 < 5/2. The validity of the predictions of Eq. ^ 
has been qualitatively verified for scale-free and scale- 
rich networks in Ref. [2l[. The physical origin of the 
different scaling for 7 larger or smaller than 5 /2 has been 
clarified in Ref. [2^]. However, a detailed investigation of 
the accuracy of the different theoretical approaches for 
generic graphs is still lacking. 

It is important to stress that Eq. © is an improve- 
ment over HMF, but it is not exact. While in HMF 
theory the actual network structure is replaced by an an- 
nealed one [23[ , QMF theory fully preserves the detailed 
quenched structure of the network as described by the 
adjacency matrix. In fact, Eq. ([T]) can be derived from 
Eq. ([2]) by inserting the largest eigenvalue of the annealed 
network, which is (k 2 )/(k) [2l[. In this sense the HMF 
approach is equivalent to QMF theory plus an additional, 
annealed network, approximation. Yet both approaches 
rely on the mean-field (single-particle) assumption that 
the probability that nearest neighbor nodes are active 
can be factorized as the product of the single node proba- 
bilities. They thus neglect possibly important dynamical 
correlations between the state of adjacent nodes. For this 
reason, while it is appropriate to say that QMF approach 
improves over HMF, there is no guarantee about the ex- 
actness of Eq. ([2]), whose accuracy needs to be checked 
on a case- by-case basis. 

In this paper we tackle this task, by performing large- 
scale numerical simulations of the SIS model on vari- 
ous types of graph, using the quasi-stationary method 
to study the phase-transition in finite systems [24J, and 
determining the effective size-dependent threshold by an- 
alyzing the peak of the susceptibility [25| . The accuracy 
of the method is checked by app lying it on heterogeneous 
annealed networks [HI, Hg, |27j . in which HMF is exact 



[23l . |28| ] . We then consider an example of homogeneous 
network, the random regular graph, and an instance of 
strongly heterogeneous network, the star network. In the 
first case, we observe that both HMF and QMF predict fi- 
nite thresholds, close to the numerical one, but not equal 
to it. In the second case, QMF correctly predicts the 
scaling of the vanishing threshold with network size, but 
the prefactor is not the right one. Finally, we consider 
the key case of power-law degree distributed networks. 
For 7 < 5/2 both theories turn out to give an asymp- 
totically exact value for the threshold. For 5/2 < 7 < 3 
we observe a vanishing threshold whose scaling with TV is 
correctly predicted by QMF, although with an incorrect 
prefactor. For 7 > 3 we numerically recover the pres- 
ence of two competing activation mechanisms discussed 
in Ref. [22j : QMF prediction follows approximately the 
transition due to the activation of the largest hub in the 
network, which vanishes and dominates in the large size 
limit. The HMF threshold remains finite and is close to 
the transition point due to the most densely connected 
core of the network (maximum fc-core). 

Our numerical results confirm that QMF is indeed an 
improvement over HMF, providing a better estimate of 
the threshold and capturing the vanishing threshold for 
power-law distributed networks with any 7, a key fact 
missed by the HMF approach. However, it is only the 
scaling of the threshold with network size that is correctly 
given by QMF. Improved analytical strategies, beyond 
quenched mean-field are thus needed in order to obtain 
more accurate threshold predictions in cases of practical 
importance. 



II. NUMERICAL METHODS 

We consider here the SIS model for epidemics in con- 
tinuous time. The numerical algorithm is implemented 
as follows: At each time step, we compute the number of 
infected nodes, TVj, and links emanating from them, N n . 
With probability Ni/(Ni + ATV„) one infected node, cho- 
sen at random, becomes healthy. With complementary 
probability ATV„/(TVj + ATV„), one of the links is selected 
and the infection is transmitted through it. The numbers 
of infected nodes and related links are updated accord- 
ingly, time is incremented by At = l/(TVj + ATV„), and 
the whole process is iterated. 



A. The quasi-stationary state method 

The standard numerical procedure to investigate the 
properties of absorbing phase transitions is based on the 
determination of the average of the order parameter (in 
this case the density of infected nodes), restricted only 
to surviving runs [8[, i.e., runs which have not reached 
the absorbing state up to a given time t. Such a tech- 
nique is quite wasteful, because surviving configurations 
are very rare at long times close to the threshold, and 
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an exceedingly large number of realizations of the pro- 
cess are needed in order to get substantial statistics. An 
alternative technique is the quasi-stationary state (QS) 
method [24], [28| , based on the idea of constraining the sys- 
tem in an active state. This procedure is implemented 
by replacing the absorbing state, every time the system 
tries to visit it, with an active configuration randomly 
taken from the history of the simulation. For this task, 
a list of M active configurations is stored and constantly 
updated. An update consists in randomly choosing a con- 
figuration in the list and replacing it by the present active 
configuration with a probability p r At. After a relaxation 
time t r , the QS quantities are determined during an av- 
eraging time t a . The QS probability P n that n vertices 
are infected is computed during the averaging interval, 
each configuration with n active vertices contributing to 
the QS distribution with a probability proportional to 
its lifespan At. From the particle distribution P n , the 
moments of the activity distribution can be computed as 
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n>l 



The simulation procedure described above was recently 
used to determine with high numerical accuracy the crit- 
ical point and exponents of the contact process [|| on an- 
nealed [28| and quenched [2!| networks with power law 
degree distributions. 

The values of the QS parameters used in the present 
simulations were M = 300 and p r = 0.02, while t r and 
t a varied depending on N and A. 



B. Numerical determination of the epidemic 
threshold 



value of A > 0, p s will attain a finite limit for a sufficiently 
large N, once the corresponding X C (N) becomes smaller 
than A. This shows that asymptotically the threshold is 



zero 



21|, but it does not provide information on the ef- 



fective threshold for finite N. To overcome this problem, 
we turn instead to another procedure to determine the 
critical point, namely the study of the susceptibility [25j . 
defined as 



(P 2 ) ~ {P? 
(P) 



(5) 



When plotted as a function of A in a system of fixed size 
the susceptibility \ exhibits a maximum at a value 
X P (N). In systems with a finite critical point in the ther- 
modynamic limit A c (oo), the peak of the susceptibility 
at X P (N) corresponds to a transition rounded by finite 
size effects, that as N — > oo tends to the critical point 
as X P (N) - A c (oo) ~ N' 1 /" @. Correspondingly, the 
height of the susceptibility at the maximum scales with 
system size as ^ max ~ A^ 7 where 7' is a new critical 
exponent. 

We assume that the relation between X p and A c written 
above holds also in the present case, where A c depends 
on N: X p (N) ~ X C (N) ~ N^ 1 ^ . This implies that the 
susceptibility peak and the size-dependent threshold tend 
to coincide in the large size limit. When the assumption 
can be explicitly controlled, it turns out (see below) to 
be correct. 

It is worth noticing that the susceptibility Eq. ([5]) is 
different from the standard definition xn = N((p 2 } — 
(p) 2 ) H- We adopt Eq. JSJ) because it leads to clearer 
numerical results (see Sec. IHIj) . while preserving all the 
scaling properties of the usual definition. 



In usual absorbing-state phase transitions, where the 
critical point converges to a finite value in the thermo- 
dynamic limit, the finite-size scaling method allows the 
precise numerical determination of the critical point and 
the whole set of associated critical exponents Q. It is 
based on the computation, for increasing system sizes N , 
of p s , the global activity density, computed from surviv- 
ing run averages or from QS simulations. For values of 
the control parameter in the active phase, p s reaches a 
finite non-zero limit as N — > 00. For values of the control 
parameter in the absorbing phase, p s decays trivially as 
p s ~ AT -1 . The critical point is distinguished by a power- 
law decay p s ~ N~P' V , where (3 is the critical exponent 
controlling the density p s at a finite distance from the 
critical point, while v is associated to the growth of cor- 
relations close to criticality ||. 

In networks with a diverging second moment for the 
degree distribution, this approach c an g o awry if strong 
corrections to scaling are present (28|.|29| . In the particu- 
lar case of the SIS model, this approach simply does not 
work, because the effective threshold depends on N and it 
goes to zero as the system size grows. Therefore, for any 



III. NUMERICAL CHECK OF THE 
SUSCEPTIBILITY METHOD: ANNEALED 
SCALE-FREE NETWORKS 

A natural benchmark to check the accuracy of the sus- 
ceptibility peak as a measure of the position of the criti- 
cal poin t in the SIS model is given by annealed networks 
[23L l26l. l27| . In annealed networks, all edges are rewired 
(preserving the degree of the involved nodes) after each 
change of the state of any vertex. This procedure de- 
stroys all types of correlations, either topological (due 
to the adjacency matrix) or dynamical, and thus renders 
exact the prediction of HMF theory [H, . From a prac- 
tical point of view, the regeneration of the whole network 
every time a microscopic dynamic step is performed can 
be effectively achieved in uncorrclated networks by se- 
lecting at random, every time that a nearest neighbor of 
a vertex is needed, a vertex of degree k', with probability 
k'P(k')/(k) [23|. In Fig. Q]we present the results of sus- 
ceptibility measurements for uncorrelated annealed net- 
works with degree exponent 7 = 2.25, and 3.5. Fig.[T](a) 
depicts the susceptibility \ defined in Eq. §5§ and the 
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FIG. 1. (Color online) (a) Plot of x = N({p 2 } - <p) 2 )/(p) 
and xn = N({p 2 } — (p) 2 ) for the SIS process on annealed SF 
networks with degree exponent 7 = 2.25 and different net- 
work sizes. The susceptibility \ is more efficient to determine 
the effective size-dependent threshold, (b) Effective threshold 
from the susceptibility peak X P (N) as a function of N, for an- 
nealed SF networks with 7 = 2.25 and 7 = 3.5. The effective 
threshold shows a very good agreement with the numerically 
evaluated threshold A™ F (Eq. ©). 



usual susceptibility \n defined in the analysis of absorb- 
ing phase transition Q, in networks with 7 = 2.25 and 
different sizes N. We observe that x provides a more 
clear-cut definition of the susceptibility peak. The same 
behavior is observed for 7 = 3.5 (data not shown). In 
Fig. Hfb) we plot the evolution of the susceptibility peak 
Ap(7V) as a function of network size for fixed 7 = 2.25 
and 7 = 3.5, and compare it with the numerically eval- 
uated HMF prediction. These results confirm that \ p 
provides an excellent approximation of the exact result 
\HMF ^ both when the threshold goes to zero with N and 
when it converges to a finite value. The differences ob- 
served might be attributed to corrections to scaling, such 
as those presented in systems with a finite threshold in 
the thermodynamic limit, see Sec. Ill Bl Therefore, in the 
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FIG. 2. (Color online) Susceptibility as function of A for RRN 
of increasing (from bottom to top) size and degree k — 10. 
The susceptibility peak is closer to the theoretical prediction 
of the pair approximation than to the HMF and QMF results. 



rest of the paper we use the position of the peak of the 
susceptibility x as the numerical estimate of the position 
of the threshold. 



IV. HOMOGENEOUS NETWORKS: THE 
RANDOM REGULAR NETWORK 



As a first non-trivial application of the technique of the 
susceptibility peak to evaluate the SIS epidemic threshold 
in homogeneous networks, we consider the case of random 
regular networks (RRN) that is, networks where all nodes 
have exactly the same degree k, while links are randomly 
distributed among them, avoiding self-connections and 
multiple connections. In this case, HMF theory predicts 
trivially a constant threshold X^ MF = 1/k. The pre- 
diction of QMF theory takes exactly the same value, as 
can be easily seen by applying Perron- Frobenius theorem 
[3fj| . Fig. [5] shows the susceptibility x as a function of A 
for RRNs with increasing N and degree k = 10. The 
numerical estimated threshold is quite off from the the- 
oretical value 1/k for HMF and QMF, indicating that 
both theories are essentially incorrect, while the suscep- 
tibility peak falls close (increasingly so for larger N) to 
the value Aj? air = l/(k — 1), which is the prediction of the 
pair approximation. 

From this analysis wc conclude that HMF and QMF 
provide a reasonable approximation but not the exact 
position of the threshold. They fail just because they 
neglect dynamical correlations among the state of neigh- 
bors, which are instead better taken into account by pair 
approximation approaches [T2I f]~3j | . 
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3. (Color online) Susceptibility \ as a function of 
ax computed in star graphs with different values of fc max 



(increasing from bottom to top) 




FIG. 4. (Color online) Effective threshold X P (N) from the 
susceptibility peak as a function of network size N for uncor- 
rected SF networks with 7 = 2.25, compared with the HMF 
and QMF predictions. Inset: Susceptibility \ as a function 
of A for different network sizes (increasing from right to left). 



V. HETEROGENEOUS NETWORKS: THE 
STAR GRAPH 

In this Section we focus on the simplest case of a het- 
erogeneous network with vanishing epidemic threshold, 
namely the star network, which is composed by a hub 
of degree fc max , to which fc max leaves of degree 1 are at- 
tached. For this star graph, the largest eigenvalue of 
the adjacency matrix can be easily shown to be vfcmax- 
The refore the QMF prediction from Eq. © is A^ MF = 
l/V^max- On the other hand, the HMF result from 
Eq. dj) takes in this case the form is A™ F = 2 /(fc ma x +l). 
Figure [3] shows the susceptibility \ versus A \J /c max com- 
puted for star graphs with a wide range of values of 
fcmax- It clearly shows that the scaling A c ~ \J k maK 
is correct, however the value of the pref actor is around 
1.5, rather than 1, in agreement with the rigorous bound 
A c > 1/V fcmax derived by Ganesh et al. [3l[. The star 
graph constitutes thus the simplest example of a net- 
work for which HMF theory does not work. This failure 
of HMF is altogether not surprising, since this particular 
network is strongly correlated at the degree level, and 
therefore fails to fulfill one necessary condition for the 
validity of the HMF result Eq. ([T]). QMF theory instead 
provides the correct scaling of the threshold with network 
size, although the prcfactor is not exact. 



VI. HETEROGENEOUS NETWORKS: 
POWER-LAW DEGREE DISTRIBUTED GRAPHS 

We now consider the SIS model on networks with 
power-law degree distributions, P(k) ~ fc~ 7 , built using 
the uncorrelated configuration model (UCM) (32|. This 
procedure is equal to the standard Molloy-Reed config- 
uration model (33l ] with the additional constraint that 
the degree values are strictly bounded by /c max ~ N 1 / 2 . 



This bound guarantees that no topological correlations 
are present in the network (34|, and therefore fulfills the 
requirement needed for the applicability of the HMF re- 
sult Eq. Q. 

We analyze three values of 7, representative of three 
regimes characterized by different expressions for the the- 
oretical estimates. 



A. 7 < 5/2 

In Fig. Uwe show the shape of the susceptibility x ver ~ 
sus A (inset) and the numerical threshold A p as a func- 
tion of the network size N (main plot) for 7 = 2.25, 
compared with the predictions of the two theoretical ap- 
proaches. It turns out that the numerical results from 
the susceptibility peak agree with good accuracy with 
both HMF and QMF theories, which in their turn tend 
to coincide. While it was expected that the theoretical 
formulas scaled in the same way with N, see Eq. ([3]), the 
fact that they tend to coincide indicates that the prcfac- 
tor c 2 in Eq. © is close to 1. Fig. [4] shows that HMF 
and QMF predictions are apparently exact in the limit 
of large systems for 7 < 5/2. 

In this particular range of 7, since the transition can 
be identified with high accuracy, it is possible to extract 
additional information about the epidemic phase transi- 
tion. We consider thus the scaling of of the quantities p s , 
X and xn with N at the transition point. According to 
standard notation Q, the expected scaling with system 
size should be (see Sec. Ill Bp : 

Ps ~N-P/\ XN~Ni'/», x „ N h'+P)/\ (6) 

In Fig. [5ja) we plot the values of p s , x an d Xn, evalu- 
ated at the susceptibility peak, as a function of N in SF 
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FIG. 5. (Color online) (a) Numerical values of p B , \ an d Xn 
as a function of TV evaluated at the susceptibility peak in SF 
networks with 7 = 2.25. (b) Stationary density of infected 
nodes as a function of the distance to the threshold A p in SF 
networks with 7 = 2.25 and TV = 10 7 . 

networks with 7 = 2.25. Fitting the curves in Fig. Eta) 
to a power law form, we find the exponents: 

/3/P = 0.65, i/9 = -0.28, (7' +P)/v = 0.37. (7) 

These exponents explain why \ is the best choice to de- 
termine the threshold. The maximum of the standard 
susceptibility \n scales with a negative exponent 7', and 
thus, in the limit of large TV, the transition is character- 
ized by a discontinuity. The value 7' + /? > instead 
ensures a clearly defined maximum for the susceptibility 
X diverging as TV — > 00. 

By plotting the order parameter p s as a function of 
the distance from the effective threshold, we can attempt 
to determine the exponent (3, which is defined by p ~ 
[A - X c (N)f. In Fig. EJb) we show such a plot, for a SF 
network with 7 = 2.25 and size TV = 10 7 . According to 
HMF theory Q, the (5 exponent is expected to take the 
value /3 = 1/(3 — 7) = 4/3, while the QMF approach of 
Van Mieghem [35| predicts /3 — 1 . The numerical results 
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FIG. 6. (Color online) Effective threshold A P (TV) as a function 
of network size TV for uncorrelated SF networks with 7 = 
2.75, compared with the HMF and QMF predictions. Inset: 
Susceptibility x as a function of A for different network sizes 
(increasing from right to left). 

presented in Fig. [SJb) yield an effective exponent lying 
between the theoretical predictions, so that the validity 
of none of them can be excluded. 



B. 5/2 < 7 < 3 

In this interval of 7, Eq. ([3]) predicts that a different 
regime sets in, with the threshold set by the inverse of 
the square-root of fc max i.e. Aj? L ~ fc ma x 2 = TV -1 / 4 , while 
HMF theory predicts A« MF ~ k m ( ^ l] = aM3-7)/2 up 
to 7 = 3. For the values of TV which can be simulated nu- 
merically, the two theoretical predictions are quite close 
but do not coincide. 

Figure[5]shows the results of the susceptibility analysis 
for SF networks with 7 = 2.75. From this plot, we con- 
clude that the numerical results do not conform to the 
HMF behavior, the more so for large system size. The 
numerical threshold A p (TV) scales instead as the inverse 
of the largest eigenvalue, but with a prefactor different 
from unity. The QMF threshold provides hence an ap- 
proximation to the numerical threshold, scaling in the 
same way, but with an accuracy of the order of 30%. 



C. 7 > 3 

For 7 > 3 HMF theory yields a finite value of the 
threshold, which instead still vanishes according to QMF. 
Since sample-to-sample fluctuations of the value of fc max 
are quite large in this regime, we consider for each value 
of TV only networks with fc max equal to the mean value 

(^max) j^lj ■ 

In Fig. [7ja) we plot the susceptibility as a function 
of A in networks with 7 = 3.5. The behavior of the 
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FIG. 7. (Color online) (a) Susceptibility \ as a function of 
A for uncorrelated networks with 7 = 3.5 and different sizes, 
(b) Effective thresholds X P (N) determined by the rightmost 
and leftmost (when present) susceptibility peak as a function 
of network size N for uncorrelated networks with 7 = 3.5, 
compared with the HMF and QMF predictions. 



susceptibility in this regime is surprisingly different from 
the one observed in the case 7 < 3. As we can see, while 
for small network sizes a well defined and unique peak 
is present for relatively large values of A, at a position 
quite compatible with the prediction of HMF, as N grows 
another feature emerges for smaller values of A, giving 
rise to a secondary peak for the largest considered sizes. 

The evidence presented in Fig. [7f a) can be better un- 
derstood with the help of the QS distribution P n of the 
order parameter (number of active nodes n = p s N) de- 
picted in Fig. [5] For large values of A, the distribution 
has a single peak (apart from the one in n = 1.). As A 
is reduced, a secondary peak starts to appear at smaller 
values of n and rapidly takes over. Further decreasing A 
leads to the disappearance of the all peaks for finite n, 
which signals the transition into the absorbing state. 

Figs. 0a) and [S] reflect the existence of two competing 
thresholds, associated to two different mechanisms trig- 
gering the transition f22j . The secondary peak for small 
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FIG. 8. (Color online) Quasi-stationary distribution of the 
number of active nodes P n for an uncorrelated network of size 
N = 3 x 10 6 and 7 = 3.5. Different curves are for decreasing 
values of A (bottom to top). 



A, whose position scales with network size as predicted by 
QMF formula, see Fig.jT^b), is associated to the presence 
of the star-subgraph centered around the largest hub, 
which for A > 1/Ajy is able to sustain alone the active 
state [2l], I22T. This kind of transition starts from a local- 
ized region [36[ and then propagates the infection to the 
rest of the network. Its position is set by fc max and does 
not change depending on the quenched network realiza- 
tion. Notice also the rounded shape of the susceptibility 
peak, reminiscent of the what is found for star graphs 
(see Fig. [3]). The peak for large A, which occurs not far 
the value predicted by HMF and is much narrower, is 
associated to the set of most densely connected nodes in 
the network (maximum k-cove) collectively turning into 
the active state. The location of this transition fluctuates 
a lot depending on the realization. 

It is clear that for asymptotically large N the first 
mechanism dominates. In this limit one expects the pic- 
ture to be analogous to the case 5/2 < 7 < 3 presented 
above: a single peak moving toward zero as N increases, 
as predicted by Eq. ^ (but with a different prcfactor). 
However, the crossover to this stage is very slow and val- 
ues of N much larger than those attainable with our com- 
putational resources would be needed to check in detail 
the accuracy of Eq. @ in this regime. 



VII. CONCLUSIONS 

In this paper we have presented a large scale numeri- 
cal analysis of the SIS model on networks. Our approach 
presents two improvements over previous numerical stud- 
ies of the SIS. Firstly, we have implemented the quasi- 
stationary state (QS) method, which allows to overcome 
the problems associated to simulations based on sur- 
viving averages, yielding far better statistics with much 
smaller uncertainties. Secondly, to overcome the prob- 
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lems that traditional finite-size scaling analysis face in 
front of a vanishing critical point, we have instead con- 
sidered the susceptibility x, whose maximum value pro- 
vides a numerical estimate X p of the threshold. The com- 
bination of the QS and susceptibility peak methods leads 
to numerically accurate threshold estimates, as we have 
checked in the case of annealed networks, in which the 
exact value of the SIS threshold is known. The accu- 
racy of our results allows us to discuss in detail the rel- 
ative performance of two candidate theoretical solutions 
to the SIS model on networks, namely the heterogeneous 
mean-field (HMF) and the quenched mean-field (QMF) 
theories. 

By considering the very simple case of the random reg- 
ular graph, our analysis shows that even for homogeneous 
networks, which have a finite threshold for large sizes, 
both HMF and QMF theories may provide inaccurate 
results. This occurs because dynamical correlations play 
a role in determining the threshold value, but they are 
disregarded by the theoretical approaches. 

In the case of strongly heterogeneous networks (the 
star graph), QMF theory is sufficient to yield the correct 
scaling, but errs in the associated prefactor, again due to 
dynamical correlations. 

Turning to the more interesting case of random un- 
corrected scale-free networks, our numerical results in- 
dicate that both HMF and QMF provide asymptotically 
exact expressions for the epidemic threshold in the case 
7 < 5/2. In the region 5/2 < 7 < 3, both theories are 
quite close for the investigated sizes but QMF is able 
to reproduce the scaling of the threshold with network 
size, erring only in the numerical prefactor. The analy- 
sis of the more complex case 7 > 3 leads to a picture in 
agreement with the presence of two epidemic activating 
mechanisms discussed in Ref . (22| . Here, the susceptibil- 
ity presents for small network sizes a peak at large A, close 
to the HMF prediction, which tends to a finite limit but 



largely fluctuates from sample to sample. This peak is 
associated to the activation of the epidemics in the net- 
work by the set of most densely connected nodes in the 
network (maximum fc-core). For large sizes, a secondary 
incipient peak appears at small A, which is described by 
QMF and asymptotically overcomes the other peak at 
sufficiently large N. The secondary peak is associated to 
the epidemic activation from the most connected node in 
the network, which as the center of a star graph of size 
+ 1 is able, all alone, to sustain activity in the whole 
network. 

The conclusion drawn from our numerical contrast of 
the HMF and QMF theories is that, while QMF repre- 
sents an notable improvement over and HMF approach, 
it is still unable to yield a precise determination of thresh- 
olds, apart from the special case 7 < 5/2. This calls thus 
for improved analytical approximations, which should in- 
clude in an explicit manner the effects of dynamical cor- 
relations between neighborin g n odes. Progress has been 
done recently along this path |l2l.[l3l.[37j. but these meth- 
ods are easily applicable only to small networks, so that 
the precise theoretical determination of the SIS epidemic 
threshold for large networks remains essentially an open 
problem. 
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